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Abstract. The critical temperature of a three-dimensional Ising model on a simple 
cubic lattice with different coupling strengths along all three spatial directions is 
calculated via the transfer matrix method and a finite size scaling for L x L x oo 
clusters (L = 2 and 3). The results obtained are compared with available calculations. 
An exact analytical solution is found for the 2 x 2 x oo Ising chain with fully anisotropic 
interactions (arbitrary J x , J y and J z ). 
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1. Introduction 

Great attention is given to the critical temperature calculation of the three-dimensional 
Ising model, for the full time extent of existence. The most considerable advances 
have been attained for the fully isotropic cubic lattice (J x = J y = J z ). The 
calculations are steadily improving with time, accuracy for the critical point value 
is down to 1(T 3 % error: K c = 0.221 6595 ± 0.000 0026 (Fcrrenbcrg and Landau 1991), 
K c = 0.221 6544 ±0.000 0010 (Livet 1991). 

The efficiency of progress is less for the partly anisotropic model when J x = 
Jy Jz- Two cases exist here: J x = J y > J z (a quasi- two-dimensional model) 
and, inversely, J x — J y < J z (a quasi-one-dimcnsional model). Using the high- 
temperature series, the critical temperature estimates have been obtained with 
2 x 10~ 2 % error for the fully isotropic interactions (J x = J y = J z ) and about 
10~ 3 % in the two-dimensional limit (J z = 0) of quasi-two-dimensional model (see 
Navarro and de Jongh 1978 and references cited there). By intermediate range of 
interlayer couplings, the error of a phase transition temperature determination lies 
between these two extreme values. Conversely, in the quasi-one-dimensional case, 
the estimates based on the same high-temperature series are rapidly deteriorated 
due to the limited number of terms available in the series. As a result, one can 
find the critical temperatures only up J x ( y )/J z = 10~ 2 (with exactness in the one- 
two significant figure range). In the quasi-one-dimensional case, the phase transition 
temperature has been calculated also by phenomenological renormalisation of clusters 
(Yurishchev and Sterlin 1991). Inasmuch, the cluster geometry reflects the physical 
situation; this approach (contrary to the high-temperature series expansions) yields 
more precise results as the anisotropy of quasi-one-dimensional system increases. By 
J x ( y )/J z = 10~ 3 , the critical temperature position is determined here with an accuracy 
of about three significant figures. 

The difficulties are largest when the interactions are different along all three 
directions. Here there are the calculations done by real-space renormalisation group 
method (da Silva et al 1984) and the calculations carried out via various versions 
of mean field theory and variational principle (see Faleiro Ferreira 1988, 1989 and 
references therein). These results we discuss in detail in the third section in comparing 
with our computations. 

In this paper, the critical temperature of a three-dimensional Ising model, with 
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fully anisotropic interactions, is calculated by the transfer matrix approach in 
combination with a finite size scaling, i.e. by the phenomenological renormalisation 
group method proposed by Nightingale (1976). As the clusters, we exploit the infinitely 
long parallelepipeds L x L x oo with transverse scales L = 2 and 3. We succeeded in 
obtaining a rigorous solution for lattice with L = 2. For the 3 x 3 x oo lattice, we 
simplify the partial eigenvalue problem for the transfer matrix of 512th order. Using 
the symmetry, we reduce this problem to a determination of the largest eigenvalues of 
the two 18th order matrices. Final calculations have already been made numerically. 



2. Calculation of the critical temperature 

Conforming with the phenomenological renormalisation group theory, the critical 
temperature T c is a fixed point of an equation (see, for instance, reviews of Nightingale 
1982 and Barber 1983): 

Lk l (T c ) = L'k l ,(T c ) (1) 

where 

k l = In(Ai/A 2 ) (2) 

is the inverse correlation length in a cluster with characteristic size L. The quantities 
Ai and A2 entering into (2) are respectively largest and next largest eigenvalues of 
subsystem transfer matrix. Thus, the task is reduced to find the dominant eigenvalues 
of transfer matrices. 



2.1. Cluster 2 x 2 x 00 

Let us write the Hamiltonian for the cluster as 

H = - y][Jx.(<Tl,j<T2,i + &3,i&A,i) + Jy{^\,i^i,i + ^2A^3,i) 
i 

+ Jz(<7l,i<7l,i+1 + &2,i&2,i+l + &3,i&3,i+l + <^i,i^i,i+l)]- (3) 

The spin variables u = ±1 arc located in sites of a lattice 2 x 2 x oo which has a 
rectangular cross section and has the symmetry planes going through its axis and the 
middles of opposite sides. 
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The transfer matrix U with elements 

(<7l, (T 2 , 0- 3 , (T 4 |J7|cri, a' 2l CT3, (T4) = CXp[iii' x ((Tl(T2 + CT 3 <74 + (T^ + CTgCT^) 

+ 5^(0-1(74 + <7 2 <7 3 + (j[(t' a + cr' 2 a' 3 ) 

+ K^a^ + a 2 a 2 + a 3 a 3 + a 4 a' i )} (4) 

where if x = J x /fcT, if y = J y /kT and iT z = J z /kT corresponds to the Hamiltonian 
(3). 

To solve the eigenvalue problem of the transfer matrix (4), we use first 
the invariance property of the appropriate Hamiltonian with respect to the 
transformations of the group Z 2 ®C2 V ; where Z 2 is a group of global reflections in the 
spin space, C2 V is the point group generated by the symmetry planes of a lattice and (x) 
represents the direct product. Carrying out the usual group-theoretical analysis (see, 
e.g., Yurishchev 1989), we come to a conclusion that the 16 x 16 transfer matrix (4) 
can be reduced due to symmetry Z. 2 ®C2 V to a quasi-diagonal form with one subblock 
5x5, four subblocks 2x2 and three 'subblocks' lxl, i.e. ready-made eigenvalues. 

The subblock of size 5 x 5 is connected with the fully symmetrical irreducible 
representation of the group. Due to the Perron theorem, it contains the largest 
eigenvalue of U. Basis vectors for this irreducible representation are given as 
. ei + eie , e 2 + e 3 + e 5 + e 8 + e 9 + e 12 + e u + e 15 

^ = ^FT ^ = 271 

(5) 

, e 4 + ei3 e 6 + en e 7 + e w 

W3 = 7= Va = 7= V5 



V2 V2 V2 



where 



ei = |1,1,1,1) e 2 = |1,1,1,-1) ... ei 6 = |-1,-1,-1,-1). (6) 

Using these basis functions and utilizing (4), we find the matrix elements ipfUtpj of 
subblock 5x5. The secular equation of this subblock has a structure (and this is the 
second key circumstance allowing the solution of the eigenvalue problem): 

A 5 - aiA 4 + a 2 A 3 - aa 2 A 2 + a 3 aiA - a 5 = 0. (7) 



Here 



ai = 2[1 + 4cosh(2if x ) cosh(2if y )] cosh(4iQ + 6 (8) 
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a 2 = 32cosh(2if x )cosh(2 J ft: y )[cosh(4is: z )cosh 2 (2/s: z ) - 1] 

+ 8[1 + cosh(4if x ) + cosh(4if y )] sinh 2 (4/C z ) (9) 

and 

a = 4sinh 2 (2^ z ). (10) 

According to Sominskii (1967), an algebraic equation like (7) is a reciprocal one. This 
property makes it possible to easily find the roots of our equation. As a result, the 
largest eigenvalue of the transfer matrix (4) is equal to 

Al = Ir 1 + (H-« 2 ) 1/2 (11) 

with 

n = i(aj - a) + [i(aj + a) 2 + a 2 - a 2 \ x l 2 . (12) 

Solving secular equations of second-order subblocks causes no difficulties. In this 
issue, we obtain a complete set of eigenvalues. Sorting the eigenvalues, we seek out 
the next largest eigenvalue of U: 

A 2 = {1 + cxp[2(if x + K y )}} sinh(4iQ + [{1 - cxp[2(K x + K y )}} 2 sinh 2 (4X z ) 

+ 16 exp[2(if x + K y )] sinh 2 (2X z )] 1 / 2 . (13) 

Note that it lies in the subblock built on basis functions which are symmetrical 
under all purely spatial transformations of the group and antisymmetrical under those 
including the spin inversion. 

By J x = J y , our solution is reduced to that of Kaufman (1949) for the Ising model 
on a cylinder, if the number of chains in the last model is equal to four. 

It is also interesting to note that the above mention does not succeed in generalizing 
the model (3). All attempts in including in the Hamiltonian new interactions (e.g., 
external field, additional pair couplings or multiparticle forces) lead immutably to the 
destruction of the obvious symmetry of Z 2 ®C2 V or the hidden algebraic one (i.e., the 
reciprocal property of a secular equation). 

2.2. Cluster 3 x 3 x oo 

We shall consider a subsystem 3 x 3 x oo with cyclic boundary conditions in both 
transverse directions. This eliminates undesirable surface effects and at the same time 
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extends the symmetry group down to Z 2 ®(T©C2 V ), where T is a group of transverse 
translations and © implies a semidirect multiplication. The given symmetry allows 
one to reduce the transfer matrix V, of the size 512 x 512, to a block diagonal form 
in which the first and second largest eigenvalues of original matrix are located in 
different subblocks (V*- 1 -* and V^ 2 \ respectively) both having a dimension of 18 x 18. 
The open form of these subblocks is given in the appendix. The extraction of dominant 
eigenvalues from and has been carried out already by the computer. 

We return again to the calculation of critical temperature. The estimates kT c /J z 
obtained by a numerical solution of transcendental equation (1) are collected in table 1. 
By this, we also put the cyclic boundary conditions on the cluster 2 x 2 x oo, i.e. simply 
increase the interaction constants in transverse directions by two times: J x — > 2 J x and 
J y — > 2J y . In table 1, we have also inserted the critical temperature values for two 
limited cases: (i) J y = 0, corresponding to the anisotropic two-dimensional Ising model 
for which the exact phase transition temperature equation is known (Onsager 1944) 



and (ii) J x = J y , corresponding to the partly anisotropic three-dimensional Ising model 
for which there exists sufficiently accurate estimates of T c (Navarro and de Jongh 1978, 
Yurishchev and Sterlin 1991). 

3. Discussion 

One of simplest ways in estimating the phase transition temperature in an Ising model 
is the mean field approximation (MFA): 



However, the accuracy is quite low (see table 2 where, for convenience of comparison, 
the critical temperature estimates found by various approximate methods have been 
given, as well as the true values obtaining from a solution of equation (14) and the 
precision numerical values). 

The state of things is somewhat corrected by an improved mean field 
approximation (IMFA), taking into account the short range order effects (Faleiro 




(14) 



(^c)mFA ~~ ^(Jx + Jy + J z )- 



(15) 
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Ferreira 1988). Inspecting table 2, one can see that there are considerable errors, 
especially for the strongly anisotropic systems. 

The MFA can be considerably improved by placing the clusters in the mean field 
instead of separate spins. A linear chain approximation (LCA) considered by Stout 
and Chisholm (1962) when the one-dimensional Ising system is taken as a cluster leads 
to the equation 

(xL— p (S) (i6 » 

where rj = (J x + J y )/J z . De Bruijn (1958) has shown the solution of an equation like 
(16) is given by 

when rj — > 0. Fisher (1967) established that the formula (17) is asymptotically exact 
for the Ising model. Although it qualitatively describes the logarithmically slow drop 
of the critical temperature with an increase of coupling anisotropy, unfortunately, 
this asymptotical formula does not provide the acceptable precision, even at high 
anisotropies. For example, by J x = J y = 1CP 2 J z a error of deviation from the high- 
temperature series estimate equals 21% and a error is 28% for the two-dimensional 
model (J y = 0) for the same value of anisotropy (J x / J z = 1CP 2 ). 

During recent years a number of equations have been obtained for the critical 
temperature of a fully anisotropic three-dimesional Ising model within the various 
generalisations of mean field theory, as well a variational approach (one from them 
- the IMFA — we have mentioned already). Using an extended variational method 
and taking the sum of linear Ising chains as an auxiliary Hamiltonian, Faleiro Ferreira 
and Silva (1982) have found an equation for T c via the so-called extended linear chain 
approximation (ELCA). The numerical solution of this equation shows that the ELCA 
perceptibly improves the LCA (see table 2), however, errors are still considerable. For 
instance, the critical point position for the two-dimensional isotropic case is overstated 
by 20%. 

Another new approach named by Faleiro Ferreira (1989), the improved linear chain 
approximation (ILCA) , is based upon the same auxiliary Hamiltonian but with other 
variational principle, leads to better results only for the isotropic three-demensional 
case (this is observed in table 2). 
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Da Silva et al (1984) have given the calculation ol a phase transition temperature 
for the fully anisotropic three-dimensional Potts lattice, a particular case of which is 
an Ising model, by the real space renormalisation group treatment. However, due to 
the finite size of spin blocks used, their estimates lose accuracy very rapidly with an 
increase of coupling anisotropy. Already only by J x = J y = 10 _1 J Z , it follows the 
overestimate kT c / J z = 1.4552 which, even after carrying out the extrapolation on the 
rather artificial scheme proposed by authors, falls to kT c /J z = 1.3986. This value 
surpasses the high-temperature series value by 4.1%. 

We now appeal to our results. The application of clusters makes it possible to take 
into consideration the specific features of a short range order and as a result, decrease 
the calculational error. Therefore, it is not surprising that the finite size scaling 
method with its hierarchy of clusters increasing in growth, allows us to determine the 
critical point of Ising model with more exactness than the approaches discussed above 
(see again a table 2). A uniform convergence (contrary to the ILCA) of the estimates 
with the growth of a lattice anisotropy is an important quality of the approximation. 
Let us consider a table 1. In the two-dimensional isotropic limit (J y = and 
Jx/Jz = 1), our calculation fixes the critical temperature with error 4.3% (in the 
direction of overestimation). For J x /J z — > 0, this error decreases continiously. This 
can be easily checked by making a comparison with the exact transition temperature 
values presented in the next column. In particular, the value has a 1.8% error by 
Jx/Jz = 10 _1 . A similar situation arises in the three-dimensional case with J x = J y . 
Here our estimates are again in excess of true values; the percentage error drops 
from the maximum value 3.8% in the fully isotropic case (J x = J y = J z ) and, for 
comparison, - to 1.7% for the J x /J z = 10 case. A analogous picture seems to 
be preserved in the intermediate region < Jy/Jx < 1: by fixing J x /J z , the error 
smoothly moves between the limited values corresponding to J y / J x = and J y /J x = 1- 

4. Conclusions 

In the present paper the more qualitative estimates of critical temperature in the fully 
anisotropic three-dimensional Ising lattice have been derived. These estimates yield 
the upper bound everywhere and by J y /J x = const the error for kT c / J z monotonically 
tends to zero when J x / J z — > 0. 
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The accurate analytical solution has been obtained for an Ising model on the 
2 x 2 x oo lattice with fully anisotropic couplings. 

The quasi-diagonalisation has been carried out for the transfer matrix of 3 x 3 x oo 
Ising model with the rectangular cross section. The expressions for the matrix elements 
of subblocks containing the leading eigenvalues are given in detail. This permits one 
to easily reproduce the results presented in the article, as well as can be useful in 
considering other problems. 
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Appendix. Explicit form of subblocks and 

The matrices and are symmetrical, therefore, it is enough to describe their 
upper triangular parts: 



V$ = 2^ \9<i' j) I cosh[( 2s - l)K z ]j exp[i(m« + m^)K x + ±( 



(Al) 



and 



y S } = 2 \/S ( T,9s'' ] sinh[(2,s - exp[|(m?+m»)/f x +i(mHm5)lf y ] (A2) 



vs=l 



where i < j = 1,2,..., 18. The basis vectors are ordered in a non-decrease of their 
lengths 

n t = {2, 6, 6, 12, 18, 18, 18, 18, 18, 36, 36, 36, 36, 36, 36, 36, 72, 72}. (A3) 

The quantities mf and m\ have a sense of the reduced partial energies of spin 
configurations in ith vector. They are 





ml 


= {9,9 


3, 3, 5, 5, 1, 1 


-3, 


1,1,-3,5,- 


3,1,-3,1 


-3} 


and 




















-{9,- 


3,9,-3,5,1,5,1 


-3, 


1,-3,1,-3, 


5,-3,1,1 


-3}. 


Finally, the weig 


at coefficients gi % '^ are g 


iven 


as 






1,1) 


1 


1,2) 


10 


2,2) 


-2 1 


1,3) 


10 


2,3) 


3 


3,3) 


-2 1 


1,4) 


10 


2,4) 


3 


3,4) 


3 


4,4) 


3 -2 1 


1,5) 


1 


2,5) 


2 10 


3,5) 


2 10 


4,5) 


4 2 


5,5) 


8 1 


1,6) 


10 


2,6) 


-2 1 


3,6) 


-1 2 


4,6) 


-2 4 


5,6) 


7 2 


6,6) 


6 2 1 


1,7) 


10 


2,7) 


-1 2 


3,7) 


-2 1 


4,7) 


-2 4 


5,7) 


7 2 


6,7) 


5 4 


7,7) 


6 2 1 


1,8) 


1 


2,8) 


2-100 


3,8) 


2-100 


4,8) 


-4 2 


5,8) 


-5 4 


6,8) 


4 -3 2 


7,8) 


4 -3 2 


8,8) 


4 -4 1 


1,9) 


1 


2,9) 


-21000 


3,9) 


-21000 


4,9) 


4 -2 


5,9) 


-5 4 


6,9) 


6-2100 


7,9) 


6-2100 


8,9) 


8 -1 


9,9) 


4 -4 1 


1,10) 10 


2,10) -1 2 


3,10) -12 



(A5) 
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4,10) -3 2 1 

8.10) 6-2100 

2.11) 1-110 

6.11) -4 4 1 

10.11) -8 8 2 

3.12) 1-110 

7.12) -4 4 1 

11.12) 10 -4 4 

3.13) -21000 

7.13) 6-2100 

11.13) -77-310 

2.14) -21000 

6.14) 6-2100 

10.14) 12 -4 2 

14.14) 3-75-21 

4.15) -23-100 

8.15) 3 -4 2 

12.15) -8 8 -2 

1.16) 1 
5,16) -5 4 

9.16) 3 -4 2 

13.16) 10 -6 2 

1.17) 10 
5,17) 6 3 

9.17) -44-100 

13.17) -77-310 

17.17) 19 -9 7 1 

4.18) -31-110 
8,18) 5-3100 

12.18) -77-310 
16,18) 9-53-10 



5,10) 7 2 

9.10) 4 -3 2 

3.11) 3 

7.11) -3 6 

11.11) 9 -5 3 1 

4.12) 2 -2 2 

8.12) -44-100 

12.12) 9 -5 3 1 

4.13) -4 2 

8.13) 4-22-10 

12.13) -12 6 

3.14) -10-110 

7.14) 2 -4 3 

11.14) -12 6 

1.15) 1 
5,15) -5 4 

9.15) 3 -4 2 

13.15) 8-44-20 

2.16) -21000 

6.16) 6-2100 

10.16) 10 -5 3 

14.16) 8-44-20 

2.17) 1-110 

6.17) -4 4 1 

10.17) -7 10 1 

14.17) -77-310 

1.18) 1 
5,18) -5 4 
9,18) 4-22-10 

13.18) 10 -6 2 
17,18) -19 13 -3 1 



6.10) 5 4 
10,10) 11 6 1 

4.11) 2 -2 2 

8.11) -44-100 

1.12) 10 
5,12) 6 3 

9.12) -52-110 

1.13) 1 
5,13) -5 4 

9.13) 5-3100 

13.13) 3-75-21 

4.14) -4 2 

8.14) 4-22-10 

12.14) -77-310 

2.15) 2-100 

6.15) 4 -3 2 

10.15) 10 -5 3 

14.15) 10 -6 2 

3.16) 2-100 

7.16) 4 -3 2 

11.16) -8 8 -2 

15.16) 11 -4 2-10 

3.17) 1-110 

7.17) -4 4 1 

11.17) 10 -4 4 

15.17) -96-210 

2.18) -21000 
6,18) 6-2100 

10.18) 8 -6 4 
14,18) 10 -6 2 
18,18) 13 -13 7 -2 1 



7.10) 5 4 

1.11) 10 
5,11) 6 3 

9.11) -52-110 

2.12) 3 

6.12) -3 6 

10.12) -8 8 2 

2.13) -10-110 

6.13) 2 -4 3 

10.13) 12 -4 2 

1.14) 1 
5,14) -5 4 

9.14) 5-3100 

13.14) 10 -6 2 

3.15) -21000 

7.15) 6-2100 

11.15) -10 4 -2 2 

15.15) 7 -8 2 1 

4.16) -23-100 

8.16) 3 -4 2 

12.16) -10 4 -2 2 

16.16) 7 -8 2 1 

4.17) 4-1100 

8.17) -52-110 

12.17) 10 -4 4 

16.17) -96-210 

3.18) -21000 
7,18) 6-2100 

11.18) -77-310 
15,18) 9-53-10 



Using (A1)-(A5) and above values for the weight coefficients, one can easily find 

12 



the expressions for the matrix elements. For instance, 
V$ = 2 cosh(9if z ) exp[9(if x + K y )] 

V { 3 \ } = 2V3[2 cosh(^ z ) + cosh(5X z )] cxp(if x + 7K y ) 

and 

Vg 18 = 2{13[sinh(if z ) - sinh(3if z )] + 7sinh(5X z ) - 2sinh(7if z ) 
+ sinh(9if z )}exp[-3(X x + K y )}. 
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Tables and table captions 

Table 1. Normalized critical temperature kT c /J z for the fully anisotropic three- 
dimensional Ising model as a function of J x /J z and J y /J x . 
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0.978 




0.892 




0.791 




0.658 


0.65} 


0.641 


0.637§ 


0.624 


0.619§ 


0.604 


0.600§ 


0.583 


0.579§ 


0.560 


0.556S 



0.004 0.431 0.428f 0.478 0.501 0.519 0.534 0.531§ 
0.003 0.410 0.407f 0.453 0.474 0.490 0.503 0.500§ 
0.002 0.383 0.380f 0.421 0.439 0.453 0.465 0.462§ 



f Onsager (1944). 

| Navarro and de Jongh (1978). 

§ Yurishchev and Sterlin (1991). 
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Table 2. Critical temperature estimates of the fully anisotropic three-dimensional 
Ising model versus calculational method. 



Met no cl 


J-X./ J~z 


= 1 


J-x./ Jz — 


io- 2 


t n 

Jy - U 


J J 


t n 

Jy - U 


J J 

<Jy J X 


MFA 


4 


6 


2.02 


2.04 


IMFA 


3.230 


4.933 


1.465 


1.487 


LCA 


3.526 


5.686 


0.590 


0.699 


ILCA 


2.885 


4.622 


0.588 


0.695 


ELCA 


2.728 


4.881 


0.543 


0.669 


Table 1 


2.367 


4.685 


0.513 


0.658 


Exact 


2.2691 


4.5115 


0.5089 


0.65 
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